## This script: individual-level correlates of local news readership using GLES 2009
## ORIGINAL DATA IS AVAILABLE FROM GESIS: https://search.gesis.org/research_data/ZA5300 
## GLES (2019). Pre-election Cross Section (GLES 2009). GESIS Data Archive, Cologne. ZA5300 Data file Version 5.0.2, https://doi.org/10.4232/1.13228.

rm(list = ls())

library(tidyverse)
library(janitor)
library(broom)
library(fixest)
library(haven)

## Custom functions 

state_id_to_names <- function(state_id){
  
  names <- recode(state_id,
                  `01` = 'Schleswig-Holstein',
                  `02` = 'Hamburg',
                  `03` = 'Niedersachsen',
                  `04` = 'Bremen',
                  `05` = 'North Rhine-Westphalia',
                  `06` = 'Hesse',
                  `07` = 'Rhineland-Palatinate',
                  `08` = 'Baden-Württemberg',
                  `09` = 'Bavaria',
                  `10` = 'Saarland',
                  `11` = 'Berlin',
                  `12` = 'Brandenburg',
                  `13` = 'Mecklenburg-Vorpommern',
                  `14` = 'Saxony',
                  `15` = 'Saxony-Anhalt',
                  `16` = 'Thuringia')
  
  return(names)
  
}

pad_zero_conditional <- function(str,
                                 str_len = 7,
                                 char_pad = '0') {
  
  ## Convert input to string, add 0, return
  str %>% as.character() %>% 
    ifelse(nchar(.) == str_len, paste0(char_pad, .), .)
  
}

## Load Data 
## Q101: Nutzung Tageszeitung: lokale oder regionale Tageszeitung

gles <- haven::read_dta('data/ZA5300_de_v5-0-2.dta') %>%
  mutate(read_ln = ifelse(q101g == 1, 1, 0),
         male = ifelse(q1 == 1, 1, 0)) %>%
  rename(age = q1a,
         pol_interest = q2,
         educ = d206) %>%
  mutate(abitur = ifelse(educ == 5, 1, 0),
         age = age/10) ## age in years divided by 10

## Code NAs

gles[gles > 91] <- NA

## Code foreign-born parents 

gles <- gles %>% 
  mutate(foreign_born_parents = ifelse(d283 == 2, 1, 0))

## Reverse code political interest

gles <- gles %>%
  mutate(pol_interest = 6 - pol_interest)

## Code Land 

gles <- gles %>%
  mutate(land = pad_zero_conditional(as.character(land), str_len = '1'),
         land_label = state_id_to_names(land),
         land_label = relevel(as.factor(land_label), ref = 'Baden-Württemberg'))

## Code foreign-born 

gles <- gles %>%
  mutate(foreign_born = ifelse(d280a > 4, 1, 0))

## Define covariates

covars <- c('male', 'age', 'abitur', 'land_label', 'foreign_born_parents', 'pol_interest')

## Create model formula 

m_formula <- paste('read_ln', '~', paste(covars, collapse = '+')) %>% as.formula()

## Estimate and tidy up

res <- fixest::feols(m_formula, data = gles) %>%
  broom::tidy(conf.int = T) %>%
  filter(!str_detect(term, 'land_label')) %>%
  filter(term != '(Intercept)') %>%
  mutate(term_label = recode(term,
                             `male` = 'Male (vs. female)',
                             `age` = 'Age (in 10-year intervals)',
                             `abitur` = 'Education: A-levels',
                             `foreign_born_parents` = 'Foreign-born parent(s)',
                             `pol_interest` = 'Political interest scale')) %>%
  mutate(term_label = fct_rev(fct_relevel(term_label, 
                                          'Male (vs. female)',
                                          'Age (in 10-year intervals)', 
                                          'Education: A-levels',
                                          'Foreign-born parent(s)',
                                          'Political interest scale')))

## Plot This 

p1 <- ggplot(res, aes(x = term_label, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_errorbar(width = 0) + 
  geom_point(shape = 21, fill = 'white') + 
  theme_bw() + 
  geom_hline(yintercept = 0, linetype = 'dotted') + 
  coord_flip() + 
  labs(x = '',
       y = 'Effect estimate')

p1

## Plot: average local news readership by state

avg_state <- gles %>%
  group_by(land_label) %>%
  summarise(read_ln = mean(read_ln, na.rm = T)) %>%
  mutate(land_label = fct_relevel(land_label, rev(levels(gles$land_label))))

p2 <- ggplot(avg_state, aes(x = land_label, y = read_ln)) + 
  geom_bar(stat = 'identity', col = 'black') + 
  theme_bw() + 
  geom_hline(yintercept = 0, linetype = 'dotted') + 
  labs(x = '',
       y = 'Share of respondents\nwho read local news')  

p2

